#################
# Appendix: USA #    
# Figure 1---3  #
#################

###load library
library(openxlsx)
library(psych)
library(ggplot2)
library(gridExtra)
library(dtw)
library(tseries)
library(TSclust)
library(tidyverse)
library(zoo)
source("extract.R")
library(fpp)



#######----------Figure 1: All patterns of indices (the US)--------##########


#######################Laakso Taagepera ENPP###########################

###-----The simple Laakso and Taagepera ENPP
usa<-read.csv("USA_FINAL.csv")
attach(usa)
usa$perdemocrat<-democrat/(democrat+republican)

#basic fomulation
usa$macropartisanship_1<-(democrat/(democrat+republican))*(effpar_ele/2)*
  (((democrat+republican)/independent)*(1/effpar_ele^2))
describe(usa$macropartisanship_1)

#weighting process
usa$macrop_var<-var(usa$macropartisanship_1)
usa$macroper_var<-var(usa$perdemocrat)

usa$macropartisanship_USA_LT<-(usa$macropartisanship_1/usa$macrop_var+usa$perdemocrat/usa$macroper_var)/
  (1/usa$macrop_var+1/usa$macroper_var)


###-------without nested s2

#basic fomulation
usa$macropartisanship_1<-(democrat/(democrat+republican))*(effpar_ele/2)*
  (((democrat+republican)/independent))
describe(usa$macropartisanship_1)

#weighting process
usa$macrop_var<-var(usa$macropartisanship_1)
usa$macroper_var<-var(usa$perdemocrat)

usa$macropartisanship_USA_LT_withoutnested<-(usa$macropartisanship_1/usa$macrop_var+usa$perdemocrat/usa$macroper_var)/
  (1/usa$macrop_var+1/usa$macroper_var)

###-----simply product the proportion of independent 


#basic fomulation
usa$macropartisanship_1<-(democrat/(democrat+republican))*(effpar_ele/2)*(independent/100)
describe(usa$macropartisanship_1)

#weighting process
usa$macrop_var<-var(usa$macropartisanship_1)
usa$macroper_var<-var(usa$perdemocrat)

usa$macropartisanship_USA_LT_prindependent<-(usa$macropartisanship_1/usa$macrop_var+usa$perdemocrat/usa$macroper_var)/
  (1/usa$macrop_var+1/usa$macroper_var)







###########################Molinar############################

###-----The simple Molinar ENPP

#basic fomulation
usa$macropartisanship_1<-(democrat/(democrat+republican))*(enpp_np/2)*(((democrat+republican)/independent)*(1/enpp_np^2))
describe(usa$macropartisanship_1)

#weighting process
usa$macrop_var<-var(usa$macropartisanship_1)
usa$macroper_var<-var(usa$perdemocrat)

usa$macropartisanship_USA_NP<-(usa$macropartisanship_1/usa$macrop_var+usa$perdemocrat/usa$macroper_var)/(1/usa$macrop_var+1/usa$macroper_var)

###-------without nested s2

#basic fomulation
usa$macropartisanship_1<-(democrat/(democrat+republican))*(enpp_np/2)*
  (((democrat+republican)/independent))
describe(usa$macropartisanship_1)

#weighting process
usa$macrop_var<-var(usa$macropartisanship_1)
usa$macroper_var<-var(usa$perdemocrat)

usa$macropartisanship_USA_NP_withoutnested<-(usa$macropartisanship_1/usa$macrop_var+usa$perdemocrat/usa$macroper_var)/
  (1/usa$macrop_var+1/usa$macroper_var)


## CV check with arima
fun2 <- function(x, h){forecast(Arima(x, order=c(1,0,0)), h=h)}
e <- tsCV(usa$macropartisanship_USA_NP_withoutnested, fun2, h=1)
sqrt(mean(e^2, na.rm=TRUE))
sqrt(mean(residuals(fun2(usa$macropartisanship_USA_NP_withoutnested,h=1))^2, na.rm=TRUE))


###-----simply product the proportion of independent 


#basic fomulation
usa$macropartisanship_1<-(democrat/(democrat+republican))*(enpp_np/2)*(independent/100)
describe(usa$macropartisanship_1)

#weighting process
usa$macrop_var<-var(usa$macropartisanship_1)
usa$macroper_var<-var(usa$perdemocrat)

usa$macropartisanship_USA_NP_prindependent<-(usa$macropartisanship_1/usa$macrop_var+usa$perdemocrat/usa$macroper_var)/
  (1/usa$macrop_var+1/usa$macroper_var)

###constituting data frame for comparison

df_mp_enpp<-data.frame(usa$yearMon,usa$perdemocrat,usa$macropartisanship_USA_LT,
                       usa$macropartisanship_USA_LT_withoutnested,usa$macropartisanship_USA_LT_prindependent,
                       usa$macropartisanship_USA_NP,usa$macropartisanship_USA_NP_withoutnested,
                       usa$macropartisanship_USA_NP_prindependent)

##################mp with DRA##############################


#load data
usae<-read.csv("usa_extract.csv")
attach(usae)
describe(usae)
usa<-na.omit(usae)

#make tidy data by tidyverse
df_usa_tidy <- usae %>%
  gather(key = VARNAME, value = value, -yearMon)
df_usa_tidy$value=df_usa_tidy$value/100

#setting time-series information
varname<-df_usa_tidy$VARNAME
date<-as.Date(df_usa_tidy$yearMon)
index<-df_usa_tidy$value
ncases<-NULL

#using "extract" to compute macropartisanship
output_usa<-extract(varname,date=date,index,ncases,begindt=ISOdate(1965,1,1),
                    enddt = ISOdate(2016,12,1), unit="M",npass=2)
display(output_usa)
summary(output_usa)

df_ld_usa<-data.frame(output_usa$varname, output_usa$loadings1,output_usa$loadings2)
colnames(df_ld_usa)<-c("Party names","Loadings1","Loadings2")
print(xtable(df_ld_usa), include.rownames=FALSE)


#factor scores of 2 dimensions
mpartisan1_USA_DRA<-output_usa$latent1
mpartisan2_USA_DRA<-output_usa$latent2
mpartisan1_USA_DRA <- mpartisan1_USA_DRA[-1]
mpartisan1_USA_DRA <- mpartisan1_USA_DRA[-623]


###constituting data frame for comparison
df_dra_usa<-data.frame(usae$yearMon,mpartisan1_USA_DRA)
colnames(df_dra_usa)<-c("usa.yearMon","mpartisan1_USA_DRA")

#join two data.frame
tb_usa_mpcomparison<-inner_join(df_dra_usa,df_mp_enpp,by=c("usa.yearMon"))
write.csv(tb_usa_mpcomparison,"tb_usa_mpcomparison.csv")


##################plot all indices#####################

#change colnames
usa_fig2<-tb_usa_mpcomparison 
colnames(usa_fig2)<-c( "usa.yearMon","DRA","DP","LT-nested","LT-simple","LT-indep","NP-nested","NT-simple","NP-independent")

#plot the simple macropartisanship and the modified
tidy_all_indices<-usa_fig2 %>% 
  gather(key="variable", value="value", -usa.yearMon)

usa_compare<-ggplot(tidy_all_indices, aes(as.yearmon(x = usa.yearMon), y = value)) + 
  geom_line(aes(color = variable), size = 1) +
  theme_minimal()
ggsave(plot=usa_compare,filename="usa_compare.pdf",width=12,height=8)



#######---------Figure 2: Dynamic time warping distance for the US case----------########

##################DTW############################

windows(30,25)
pdf("fig_dtw_usa.pdf")
split.screen(c(3,3))


###----displaying dynamic time warping:

####DRA
#simply calulating dtw distance
df_dtw<-data.frame(tb_usa_mpcomparison$usa.perdemocrat ,tb_usa_mpcomparison$mpartisan1_USA_DRA)
df_dtw<-na.omit(df_dtw)

ts_a<-df_dtw$tb_usa_mpcomparison.usa.perdemocrat
ts_b<-df_dtw$tb_usa_mpcomparison.mpartisan1_USA_DRA

d2 <- dtw::dtw(ts_a, ts_b,step.pattern = symmetric1)
d2$distance

windows(10,8)
dt7<-dtwPlotTwoWay(d2, ts_a, ts_b,main="CONVENTIONAL vs DRA")


###LT nested
#simply calulating dtw distance
df_dtw<-data.frame(tb_usa_mpcomparison$usa.perdemocrat ,tb_usa_mpcomparison$usa.macropartisanship_USA_LT)
df_dtw<-na.omit(df_dtw)

ts_a<-df_dtw$tb_usa_mpcomparison.usa.perdemocrat
ts_b<-df_dtw$tb_usa_mpcomparison.usa.macropartisanship_USA_LT

d2 <- dtw::dtw(ts_a, ts_b,step.pattern = symmetric1)
d2$distance

windows(10,8)
dt3<-dtwPlotTwoWay(d2, ts_a, ts_b,main="CONVENTIONAL vs LT-nested")


###NP nested
#simply calulating dtw distance
df_dtw<-data.frame(tb_usa_mpcomparison$usa.perdemocrat ,tb_usa_mpcomparison$usa.macropartisanship_USA_NP)
df_dtw<-na.omit(df_dtw)

ts_a<-df_dtw$tb_usa_mpcomparison.usa.perdemocrat
ts_b<-df_dtw$tb_usa_mpcomparison.usa.macropartisanship_USA_NP

d2 <- dtw::dtw(ts_a, ts_b,step.pattern = symmetric1)
d2$distance

windows(10,8)
dt4<-dtwPlotTwoWay(d2, ts_a, ts_b,main="CONVENTIONAL vs NP-nested")


###LT simple
#simply calulating dtw distance
df_dtw<-data.frame(tb_usa_mpcomparison$usa.perdemocrat ,tb_usa_mpcomparison$usa.macropartisanship_USA_LT_withoutnested)
df_dtw<-na.omit(df_dtw)

ts_a<-df_dtw$tb_usa_mpcomparison.usa.perdemocrat
ts_b<-df_dtw$tb_usa_mpcomparison.usa.macropartisanship_USA_LT_withoutnested

d2 <- dtw::dtw(ts_a, ts_b,step.pattern = symmetric1)
d2$distance

windows(10,8)
dt1<-dtwPlotTwoWay(d2, ts_a, ts_b,main="CONVENTIONAL vs LT-simple")


###NP simple
#simply calulating dtw distance
df_dtw<-data.frame(tb_usa_mpcomparison$usa.perdemocrat ,tb_usa_mpcomparison$usa.macropartisanship_USA_NP_withoutnested)
df_dtw<-na.omit(df_dtw)

ts_a<-df_dtw$tb_usa_mpcomparison.usa.perdemocrat
ts_b<-df_dtw$tb_usa_mpcomparison.usa.macropartisanship_USA_NP_withoutnested

d2 <- dtw::dtw(ts_a, ts_b,step.pattern = symmetric1)
d2$distance

windows(10,8)
dt2<-dtwPlotTwoWay(d2, ts_a, ts_b,main="CONVENTIONAL vs NP-simple")

###LT indep
#simply calulating dtw distance
df_dtw<-data.frame(tb_usa_mpcomparison$usa.perdemocrat ,tb_usa_mpcomparison$usa.macropartisanship_USA_LT_prindependent)
df_dtw<-na.omit(df_dtw)

ts_a<-df_dtw$tb_usa_mpcomparison.usa.perdemocrat
ts_b<-df_dtw$tb_usa_mpcomparison.usa.macropartisanship_USA_LT_prindependent

d2 <- dtw::dtw(ts_a, ts_b,step.pattern = symmetric1)
d2$distance

windows(10,8)
dt5<-dtwPlotTwoWay(d2, ts_a, ts_b,main="CONVENTIONAL vs LT-indep")


###NP indep
#simply calulating dtw distance
df_dtw<-data.frame(tb_usa_mpcomparison$usa.perdemocrat ,tb_usa_mpcomparison$usa.macropartisanship_USA_NP_prindependent)
df_dtw<-na.omit(df_dtw)

ts_a<-df_dtw$tb_usa_mpcomparison.usa.perdemocrat
ts_b<-df_dtw$tb_usa_mpcomparison.usa.macropartisanship_USA_NP_prindependent

d2 <- dtw::dtw(ts_a, ts_b,step.pattern = symmetric1)
d2$distance

windows(10,8)
dt6<-dtwPlotTwoWay(d2, ts_a, ts_b,main="CONVENTIONAL vs NP-indep")



#######---------Figure 3 Dynamic time warping distance for the US case----------########

##################TSCA############################

#change colnames
summary(tb_usa_mpcomparison)
df_tsclust<-tb_usa_mpcomparison[,-1]
colnames(df_tsclust)=c("DRA","CONVENTIONAL","LT-nested","LT-simple","LT-withIndep",
                       "NP-nested","NP-simple","NP-withIndep")


# calculating dtwarp distance
d_dtw_clust <- diss(df_tsclust, "DTWARP")
d_dtw_clust
# hclust with 最遠隣法
h <- hclust(d_dtw_clust)

windows(10,8)
pdf("fig_tclust_usa.pdf")
par(cex=0.8)
plot(h,main="Time-series Cluster Dendrogram" ,xlab="cluster(complete method)",hang = -1)
dev.off()



